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Abstract. Turbulence in compressible plasma plays a key role in many areas of 
astrophysics and engineering. The extreme plasma parameters in these environments, 
e.g. high Reynolds numbers, supersonic and super-Alfvenic flows, however, make direct 
numerical simulations computationally intractable even for the simplest treatment - 
magnetohydrodynamics (MHD). To overcome this problem one can use subgrid-scale 
(SGS) closures - models for the influence of unresolved, subgrid-scales on the resolved 
ones. In this work we propose and validate a set of constant coefficient closures for the 
resolved, compressible, ideal MHD equations. The subgrid-scale energies are modeled 
by Smagorinsky-like equilibrium closures. The turbulent stresses and the electromotive 
force (EMF) are described by expressions that are nonlinear in terms of large scale 
velocity and magnetic field gradients. To verify the closures we conduct a priori tests 
over 137 simulation snapshots from two different codes with varying ratios of thermal 
to magnetic pressure (/3 P = 0.25,1, 2.5, 5, 25) and sonic Mach numbers (M s = 2, 2.5,4). 
Furthermore, we make a comparison to traditional, phenomenological eddy-viscosity 
and a — /? — 7 closures. We find only mediocre performance of the kinetic eddy- 
viscosity and a — /3 — 7 closures, and that the magnetic eddy-viscosity closure is 
poorly correlated with the simulation data. Moreover, three of five coefficients of the 
traditional closures exhibit a significant spread in values. In contrast, our new closures 
demonstrate consistently high correlation and constant coefficient values over time and 
and over the wide range of parameters tested. Important aspects in compressible MHD 
turbulence such as the bi-directional energy cascade, turbulent magnetic pressure and 
proper alignment of the EMF are well described by our new closures. 
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Turbulence is ubiquitous in astrophysical plasmas, ranging from coronal mass ejections 
and stellar winds [1], through star formation in molecular clouds [2], to the gas in 
the interstellar [3] and intracluster medium. While experimental setups [4] become 
increasingly more realistic, they are still far away from the regime acting in such 
extreme conditions. For numerical simulations it is computationally too expensive (if 
even possible) to capture the entire range of physical processes from plasma kinetics 
to the integral scales of turbulence. In an astrophysical context, one has to further 
contend with the additional complications brought about by high compressibility and 
the accompanying supersonic and super-Alfvcnic motion. 

Possible ways to circumvent the infeasibility of direct numerical simulations are the 
use of calculations based on mean-field theories or large-eddy simulations [5]. These 
simulations only resolve the energy containing large scale dynamics and require a 
subgrid-scale (SGS) model to account for unresolved effects. While a lot of research 
has been successfully carried out in the realm of hydrodynamics [6], compressible 
magnetohydrodynamic SGS closures are essentially unexplored. Previous research is 
mainly based on the concept of turbulent dissipation in incompressible flows [7, 8, 9, 10]. 
They expand the idea of a turbulent eddy-viscosity to an additional eddy-resistivity in 
the induction equation and propose different phenomenological models. Even though 
these models are then evaluated a posteriori , a general verification and justification 
a priori has so far only been considered for a single incompressible dataset [11]. 
Thus, our objective is to establish the validity of closures for the filtered, compressible 
MHD equations by coarse-graining multiple datasets from high-resolution simulations 
of statistically homogeneous, forced MHD turbulence. 

In general, the effect of finite resolution in numerical simulations can be mimicked 
by applying a low-pass filter to the standard, ideal MHD equations. This is achieved by 
convolving the equations with a suitable filter kernel G. See e.g. Gamier et al [6] for 
details on the properties of low-pass filtering and the conditions that G needs to satisfy. 
For a homogeneous, isotropic, stationary kernel, under periodic boundary conditions 
[12] the equations take the following form 

! + v.(pa) = o, 


^—b V ■ ( pu — B ® B ) + V 

at v 7 

d~B 

— -V x (u x B) = V x £. 

dt K ’ 



-V-r, 


The units of the magnetic held B incorporate l/-\/47i\ An overbar □ denotes^ filtered 
and a tilde □ mass-weighted filtered quantities [12]. For instance, the filtered density 


§ Throughout the paper the symbol □ is used as a generic placeholder for variables. □ designates 
closure expressions. Furthermore, we employ Einstein summation convention and is identified 
with the fc-th partial derivative of the *-th component of □. 
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field is given by p = G* p, while the mass-weighed filtered velocity held is u = up/~p. In 
this formalism, all filtered primary quantities, density p, velocity u, magnetic held B , 
and thermal pressure P are presumed to be known and directly accessible. Due to the 
introduction of mass-weighted hltering the only remaining terms that require closure are 
the SGS stress r and the electromotive force (EMF), S. They are analytically expressed 
[10] as 

£ = u x B — u x B and (1) 

pj = r ij - r g + ( B2 - ^ 2 ) y, with 

rfj = ~p(upfj - UiUj) and T b - = (BiBj - B, B 3 ) . (2) 

The SGS stress tensor can be decomposed into the well-known turbulent Reynolds stress 
r u , a turbulent Maxwell stress r b and a magnetic pressure term. 

Furthermore, the definitions of the SGS energies are obtained from applying 
the hlter to the total hltered energy density E, which can be decomposed into the 
contribution due to resolved helds only and a remainder, designated as SGS energy 

E=^ + ^ 2 + lK“ 2 -“' 2 ) + i(^- 52 )' 

S - V -' V -V-' 

(resolved) =E“ ga +E| > gg =E aga (unresolved) 

It is important to point out that in general the hltering operator is not a Reynolds 
operator, in particular □ ^ It follows that SGS terms, like E$gs , carry information 
not only about the interactions between unresolved helds but also about cross-scale 
interactions between unresolved and resolved helds. In addition to this, the turbulent 
magnetic pressure is identical to the magnetic SGS energy F b gs and both kinetic and 
magnetic SGS energies are directly given by 2Ef gs = Tr (r u ) and 2 E^ = Tr (r b ), i.e. 
they constitute the isotropic parts of the respective SGS tensors. Following the general 
tensor decomposition, the deviatoric, traceless parts are then given by r,fj* = Rj - 1 ^ij T kk- 

2. Traditional closures 

In hydrodynamics, the traceless part of the SGS stress tensor is commonly closed 
by means of the eddy-viscosity hypothesis rf* = —2u u pS* 3 in analogy to the 
molecular viscosity term in the momentum equation, where S 7J = \ (u l:J + Uj t f) is the 
hltered kinetic rate-of-strain tensor. This introduces a turbulent kinetic eddy-viscosity 
u n = Cf A (E” /p) ” which is proportional to a characteristic velocity, commonly 
given by the kinetic SGS energy, and a characteristic length scale A. This closure has 
already been applied directly to MHD [7, 13] by neglecting the magnetic contribution 
rf* in the momentum equation. A turbulent magnetic viscosity u h = C* b A (F b gs ) 12 was 
used in [10] with the closure 

rg* = -2 u h M~ 3 , 

where N4 tJ = | (R, ;j + B 3 f) is the hltered magnetic rate-of-strain tensor. 
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The SGS energies can either be determined by individual evolution equations, where 
several terms again require closure, or by an instantaneous closure. Smagorinsky [14] 
introduced such an instantaneous closure in pure incompressible hydrodynamics ( B = 0) 
by assuming the SGS energy flux to be in equilibrium with the rate of dissipation 

E^ = C£A 2 p\S*\ 2 . (3) 

Here, |«S*| = \J 2 <SfSf denotes the rate-of-strain magnitude. 

Finally, the EMF is commonly modeled (e.g. [7, 8, 13, 10]) by variations of [15] 

£ = aB — (3J + yfi, 

with resolved current J = V x B and vorticity f2 = V x u. The coefficients a, (3, and 
7 are typically related to the a-effect, turbulent resistivity and turbulent cross hclicity, 
respectively. The commonly used closures for these coefficients, 

& G^tturb-Tf, f G^jtturb-S'sgs /Pi '7 C^fturbkF, 

are based on dimensional arguments, with turbulent cross hclicity W = u • B — 
u ■ B, residual hclicity H ~ (j • B-JB) — p • u — • uj , and timescale 

tturb = A (E sgs /p)~ 1/2 . 


3. Nonlinear closures 


In our new approach we adopt the compressible hydrodynamic nonlinear closure for 
the kinematic deviatoric stress tensor r u * from [16], similar to the incompressible one 
from [17]. We propose the straightforward extension to MHD with 


ij 


Ti 


b* 


or 1 ' 1 j? u 

^'nWsgs 


o/^b 7-ib 

ZL nl- C/ sgs 


1 \ 

\ Ui^s 3 J 
\B liS B l)S 3 dlJ ) 


(4) 

(5) 


The tensorial structure, e.g. hyfcUjy, can be obtained by a Taylor expansion 
discarding terms with 2 nd and higher order gradients of the resolved fields. The overall 
normalisation with the subgrid-scale energies comes from the constraint that the SGS 
stresses vanish in laminar flows. 

Applying the nonlinearity idea to the EMF generalizes the closure proposed by [11] 
to the compressible regime 

&ijkC jjA Uj tS Bk,s ■ (6) 


The closure explicitly preserves the anti-symmetry between velocity and magnetic field 
in £, which in turn helps in capturing their relative geometry. 

Finally, to complete the set of nonlinear closure equations, we use the Smagorinsky 
expression for the turbulent kinetic energy (3) and propose an analogous extension to 
the magnetic part 

E s b gs = G b A 2 |AT| 2 • 


(7) 
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Here, the turbulent magnetic energy is proportional to the magnetic rate-of-strain 

magnitude \M.\ = \J2JA l3 A4,j. There are two advantages of closing Ef gs and E b gs 
separately and not jointly via the total SGS energy. First, there is no additional need to 
close the often neglected turbulent magnetic pressure, as it is given by E b gs . Second, the 
individual energies provide closures to the isotropic parts of the turbulent stress tensors 
r u and r b . 

4. Validation Method 

In order to evaluate the proposed closures, we perform an a priori comparison 
using simulation data obtained from two, grid-based MHD codes (Enzo [18] and 
FLASHv4 [19]). This way the results are less likely to hinge on the particulars of 
the numerical implementation. In both cases we follow the evolution of a compressible, 
isothermal fluid in a cubic box with resolution of 512 3 grid cells and periodic boundary 
conditions, starting from uniform initial conditions. In Enzo we use an ideal equation 
of state with adiabatic exponent k = 1.001 in order to approximate isothermal gas. 
Enzo is a finite-volume code, i.e. the evolution equations are evaluated in integral form 
by solving a Riemann problem for the mass, momentum and energy flux through cell 
walls. This allows for the conservation of MHD invariants (e.g. energy) to machine 
precision. We use a MUSCL-Hancock scheme [20] (a second-order accurate Godunov 
extension) with second-order Runge-Kutta time integration and Harten-Lax-van Leer 
(HLL) Riemann solver (a two-wave, three-state solver) to solve the ideal MHD equations. 
The FLASHv4 code is similar to Enzo (second-order accurate in space and time), 
but uses the positive-definite HLL3R Riemann solver [20]. Another difference is that 
FLASHv4 uses a polytropic equation of state to keep the gas exactly isothermal. 
Moreover, explicit kinematic viscosity and magnetic resistivity terms are included in 
the momentum, energy, and induction equations. We set the kinematic and magnetic 
Reynolds numbers to Re = Rm = 3780. Consequently, the magnetic Prandtl number 
Pm = Rm/Re is unity. For details on the numerical methods used in FLASHvT, 
including viscous and resistive dissipation, see [21] and [22], Both codes employ 
divergence cleaning [23] to maintain V • B — 0. A state of homogeneous and isotropic 
turbulence is reached by supersonic stochastic driving in the momentum equation (given 
by an Ornstein-Uhlenbeck process) at small wave-numbers, similar to [24, 25]. Thus, 
the forcing field is evolving in time and space. The associated large auto-correlation 
time-scale T of the forcing translates to the eddy turnover time of the largest, energy- 
containing eddies. It is therefore the chosen unit of time in the following. 

We explore a range of parameters. The initial strength of the magnetic field is set 
by the plasma /3 P - the ratio of thermal to magnetic pressure. The final sonic Mach 
number M s is determined by the forcing amplitude. For the Enzo simulations we have 
initial /3 P = 0.25, 2.5, 25, with M s 2.5 after t m 2T turnover times. The FLASHv4 
simulations reach M s « 4, 2 for initial /3 P = 1,5, keeping instead constant Alfvenic Mach 
number M a 3. We discard all initial data affected by transients (before a simulation 
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time of t — 2 T) and analyze consequent snapshots taken approximately in intervals of 
0.15T and 0.1T for the Enzo and FLASHv4 datasets, respectively. 

The analysis begins with the application of a low-pass Gaussian (test) filter to the 
equations of motion. In the context of the closures we investigate, filtered quantities (i.e. 
density, velocity, and magnetic held) are interpreted as resolved, while the remainders 
represent the unresolved small scales. We can then compute r and £ both directly from 
(1) and (2), and from their respective traditional and nonlinear closures. 

The determination of the length-scale of the filter bears some consideration. It 
needs to fall within an intermediate range of length scales, away from the particular 
effects of both the large-scale forcing and the small-scale dissipation. The largest scale 
of the system is the full box size L and corresponds to the wavenumber n — 1, while 
the smallest scale is given by the Nyquist wavenumber %y q = N/2 = 256 for the 
linear numerical resolution N = 512 grid cells. The turbulence injection wavenumber is 
ninj = 2 (corresponding to half the box size L/ 2) in both codes, which is why the energy 
spectra, as illustrated in figure 1, peak there. The figure shows the mean kinetic and 
total (kinetic plus magnetic) energy spectra as a function of wavenumber. Since the 
stochastic forcing is implemented only in the momentum equation both for Enzo and 
FLASHv4, the kinetic energy spectrum exhibits the most direct imprint of the forcing 
itself. Conversely, the total energy spectrum carries the overall effect of the small- 
scale dissipation through both kinetic and magnetic channels. Figure 1 demonstrates 
that our simulations produce approximate power-law scaling within a narrow range of 
wavenumbers, which is indicative of self-similar turbulent fluctuations [26]. This can 
be interpreted as inertial range dynamics, although the nature of the inertial range in 
compressible MHD turbulence is still not fully understood [27, 28, 29, 30]. Furthermore, 
this range separates the forcing scale and the dissipation scales and is not affected by 
numerical diffusion in the absence of a bottleneck effect as demonstrated by [31]. The 
vertical dotted line in figure 1 indicates our chosen filter length scale, corresponding to 
A = 16 grid cells or wavenumber nfii ter = N/ (2A) = 16. This filter scale falls within the 
range of the self-similar power-law range for both the kinetic and total energy spectra. 
This is why we use this ideal scale for our filter in the following analysis. 

Additionally, this provides the motivation to treat data from both simulations on 
equal footing, even though FLASHv4 has explicit viscosity and diffusivity while Enzo 
solves the ideal MHD equation, subject to numerical dissipation only. 

In order to incorporate coordinate independence, a scalar held is chosen for 
comparison, the SGS energy flux E, i.e. the term responsible for the transfer of SGS 
energy between resolved and unresolved scales. Its components associated with the 
Reynolds and Maxwell stresses and the EMF are = rfjSij, and = £■ J, 

respectively (see appendix of [10] for the detailed SGS energy equation). Here, we 
substitute (1) and (2) to obtain exact fluxes E D and match these to model fluxes E n 
that employ the corresponding closures. For example, in the case of the eddy-viscosity 
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(a) kinetic energy 


(b) kinetic + magnetic energy 


Figure 1. Kinetic (a) and total (b) energy spectrum for each dataset, averaged over the 
time between 2 T and 5 T. The kinetic energy is calculated from the Fourier transform 
of yfpu. 


closure for the deviatoric turbulent Reynolds stress tensor we compare the exact flux 

S Tjj Sij p (^UiUj 'Uj'itj) Sij 

with the model flux 


( 8 ) 


£"• = = -Vp|S'| 2 = -c ;A (pE“P ) 1/2 |S-| 2 . (9) 

On the one hand, the comparison involves the determination of the constant (in space 
and time), dimensionless closure coefficients Cff. They are computed individually for 
each snapshot by minimizing the error between S D and X D in the least-square sense. 
This allows to further test the constancy of the coefficients with respect to time and 
plasma parameters. On the other hand, the general performance of the closure is gauged 
by computing the Pearson correlation coefficient of S D and £ n , where the obtained 
closure coefficients are substituted in. 

Several assumptions should be pointed out concerning this validation technique. 
Firstly, the simulation data we have available for comparison fall short of realistic 
astrophysical parameters, e.g. with regards to Reynolds numbers and resolution. In 
that sense, it would be interesting to use higher resolution direct numerical simulation 
data or three-dimensional observations or experimental results. The problem is that 
experimental data for supersonic compressible turbulent plasmas are not available and 
obtaining realistic Reynolds numbers is computationally challenging for astrophysical 
parameters. However, as seen from figure 1, the data we have are sufficiently well 
resolved for our analysis. Secondly, in choosing the SGS energy flux £, as a diagnostic 
variable, we implicitly assume that in the context of homogeneous and isotropic 
turbulence the turbulent transport (encoded by terms of the form V • {u ■ r) and 
V ■ (yB x £)) averages out to zero on subgrid scales. This assumption can nevertheless 
be easily relaxed by incorporating further diagnostic variables. Finally, we have focused 
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Table 1 . Model coefficient overview - coefficient value and energy flux correlation: 
median and bounds of the central 90% interval across all datasets. 


Model 

Coefficient 

Value 

Corr[E n , E n ] 

Smagorinsky 

G£ 

0-056t°;°l 6 5 

0.075™ 

n qn+ 0 . 02 4 
u,b,u -0.04 

0 qi +0.021 

u - yi -0.04 

eddy-viscosity 

Cf 

G b 

0 O61 +0 ' 045 

u.uui_ 0 019 
- 0 - 002™ 9 

0 70+ 013 
u - /u - 0.11 

0 06 +o ' 14 
u,uu -0.06 

nonlinear 

G n u , 

G n b i 

Gf, 

0 . 68 ™ 
n 77+O.O8 

U. / l —0.12 
n 1 O+0.013 
u * lz -0.024 

n 94+ 0 - 04 

u - y4: -0.04 

0 QQ+0.04 

u -^ u -0.07 

0.7918;?? 

a — /3 — 7 

Gf 

Gf 

Gf 

0 0007 +0 ’ 0010 

U.UUUf _ 0 0016 

0 020 + °' 009 
U.UZU _ 0 005 

—0 OO ^ 0,067 
U.UUO_o 045 

| 0.58™ 


on the SGS energy since it increases monotonically with the strength of turbulence 
regardless of the type of turbulence (e.g. compressive or solenoidal, weak or strong, etc.). 
As an extension, the other two quadratic MHD invariants - the magnetic helicity and 
cross-helicity, may further highlight distinct turbulence properties present in particular 
flow configurations. These should be kept in mind as further avenues of investigation, 
once a preferred closure has been identified by the described validation technique. 

5. Results 

The fitting results for the SGS stress tensors’ energy flux are given in figure 2 for the 
isotropic components and figure 3 for the deviatoric components. 

The isotropic parts of r u and r b are given by the SGS energies r? = |£P s from 
(3) and (7). Both, the coefficient values of kinetic part (figure 2(a) top panel) and 
the magnetic part (figure 2(b) top panel), have a small spread within a factor of two 
across time and all simulations. Furthermore, closure and data are highly correlated 
(bottom panels) with a median correlation coefficient of 0.90 and 0.91, respectively. 
More detailed numerical values of these and all following coefficients and correlations 
are listed in table 1. 

The differences in the deviatoric parts r u * in figure 3(a) and r b * in figure 3(b) 
between the nonlinear and the eddy-viscosity closures are apparent. While our nonlinear 
closure exhibits approximately constant coefficient values and correlations over time in 
all simulations, the kinetic eddy-viscosity closure shows a correlation weaker by ~ 0.2 
and bigger spread in coefficient values. Moreover, the magnetic eddy-viscosity closure 
is effectively uncorrelated with the simulation data and the coefficients can even switch 
sign at different times. The performance of the different closures can be understood from 
figure 4, where we plot the energy flux distributions E u * and E b * for a single snapshot. 
A negative flux E u * < 0 corresponds to a forward energy cascade — the transfer of energy 








Corr £ 
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1 L5 
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=Th 0.9 

sh 0.6 

, M i 

a o.3 

o 

O 


Smag. closure 

° kin. 

i 

2.5 3.5 4.5 

time t/T 




i 

2.5 3.5 4.5 

time t/T 


(a) kinetic SGS energy 


(b) magnetic SGS energy 


Figure 2. Model coefficient values (top panels) and correlations (bottom panels with 
— median) from fitting model energy flux Eg = rf/Stj to exact flux for the isotropic 
parts of the SGS stress tensors. These are given by the respective energy model in the 
trace elements (tP = | /y^ s ). Each panel contains the joint data of all simulations and 
each snapshot is represented by a marker. Values are given in table 1. 




(a) traceless kinetic SGS stress tensor (b) traceless magnetic SGS stress tensor 

Figure 3. Model coefficient values (top panels) and correlations (bottom panels with 
median) from fitting model energy flux E n = rf/Sij to exact flux for the nonlinear 
closure (left panels) and eddy-viscosity closure (right panels). Each panel contains the 
joint data of all simulations and each snapshot is represented by a marker. Values are 
given in table 1. 
























































Nonlinear closures for scale separation in supersonic MHD turbulence 


10 



Figure 4. Representative snapshot (Enzo sim. with /3 P = 2.5 at t = 4.44T) of the 
energy flux = T ij^ij distribution within the simulation box. 

from resolved to subgrid scales, because £ u * appears as a sink term in the SGS kinetic 
and magnetic energy evolution equations and as a source term in the respective resolved 
energy equations. Conversely, a positive flux corresponds to an inverse energy cascade, 
i.e. transport of energy from subgrid to resolved spatial scales. The general distribution 
of the actual fluxes in figure 4 is representative for all snapshots. The kinetic SGS energy 
fluxes are globally almost 1:1 in both directions of the turbulent cascade with a slight 
tendency towards the forward cascade. However, the forward cascade is about 3-10 
times stronger depending on the parameters as indicated by the position of the peaks 
in the distribution. For this reason, the kinetic eddy-viscosity closure shows a moderate 
correlation even though it captures only the forward energy cascade - from large to 
small scales. In fact, since under the eddy-viscosity hypothesis the kinetic SGS energy 
flux has the form £ u * = — u u J)\S*\ 2 , see (9), any model in which the eddy-viscosity u u 
has a de fini te signature with respect to space cannot reproduce a bi-directional energy 
cascade that is well represented by the nonlinear closure. In contrast to the kinetic 
SGS energy flux, the global magnetic flux clearly has a preferred direction. Depending 
on the parameters, between 60% and 80% of cells have a positive SGS magnetic flux 
indicating energy transfer from large to small scales. Nevertheless, the difference in 
strength is less pronounced as the overall forward flux is about 2 times stronger than 
the backward one. Again, these properties are well captured by the proposed nonlinear 
closure whereas the magnetic eddy-viscosity closure is poorly correlated in both strength 
and magnitude. Finally, it should be noted, that the nonlinear closures also work very 
well with the Smagorinsky energy closure. Exchanging the exact expressions E^ gs in (4) 
and (5) with £% s only slightly reduces the correlations (max. 5%) and the coefficients 
remain constant up to the second significant figure (not plotted here). 

Moving on to the EMF, the nonlinear closure outperforms the traditional ct — /3 — 7 
closure in almost all datasets, maintaining a constant coefficient with a median 
correlation of 0.79 (figure 5). The traditional closure exhibits consistently weaker 
correlations, despite the increased flexibility of three free coefficients. Only Cf, related 
to the turbulent resistivity term in the EMF, is approximately constant, whereas Cf and 
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time t/T time t/T 


Figure 5. Model coefficient values (top panels) normalized to the sample median 
(—) and the corresponding Pearson correlation coefficients (bottom panels) with 90% 
central interval (- -) for the nonlinear closure (left panels) and the reference closure 
(right panels) from energy flux fitting for £. Each panel contains the joint data of all 
simulations and each snapshot is represented by a marker. Values are listed in table 1. 


C!r fail to maintain steady values or consistent signature. The reason for the consistently 
better correlations of the nonlinear closure is hinted at in figure 6. This probability 
density plot of the local alignment between S and S demonstrates that the traditional 
closure is almost randomly aligned (flat distribution) whereas the nonlinear closure 
approaches the desired ^-distribution at 0°. 

6. Conclusions and outlook 

In summary, we have proposed a set of constant coefficient closures for the SGS 
stress and EMF in the filtered MHD equations and conducted a priori tests. The 
tests we performed do show that the new nonlinear closures perform significantly 
better than traditional, phenomenological closures with respect to both structural and 
functional diagnostics. The tests consist of filtering Enzo and FLASHvd simulations 
of homogeneous, isotropic turbulence and comparing the resulting SGS terms to their 
respective closures (dependent only on the hltered fields). All quantities are compared 
via their contributions to the SGS energy flux E n , where the closure coefficients 
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150 120° 90° 60 300° 



Figure 6. Probability density plot of the local (cell-by-cell) alignment between £ and 
£. The lines designate the median PDF across all times and datasets for the nonlinear 
closure (red —) and the a — j3 — 7 closure (gray - -). The shaded regions designate the 
central 90% interval across all times and datasets. 

are computed by individual least-square fitting. In addition, the alignment for the 
EMF vector is investigated. All new coefficients correlate well with the data. They 
are constant over time and as a direct consequence the proposed closures may be 
implemented in large-eddy simulations without the need for a computationally expensive 
dynamical procedure which computes the coefficient values at run time. In addition, the 
coefficients remain constant across simulation runs from two different codes and a wide 
range of plasma parameters, suggesting that the proposed closures capture an underlying 
physical mechanism at work in highly compressible turbulent plasma flows. Moreover, 
the new closures successfully represent the turbulent magnetic pressure, reproduce the 
bi-directional energy cascade and are well aligned with the EMF. We recognize the 
slightly lower correlation of the nonlinear closure in the EMF than in the SGS stress 
counterpart, suggesting small room for improvement. 

Nevertheless, the performance improvement over the traditional closures already 
supports the implementation and validation of the new closures in an SGS model for 
large-eddy simulations of compressible turbulent plasma flows. These simulations would 
then allow us to infer the effect of the proposed model on the large scale flow in 
practice. Potential applications include accretion disks [32], star-forming magnetized 
clouds [33, 34] and plasmas on cosmological scales [35, 36, 37, 38, 39, 40]. 



























Nonlinear closures for scale separation in supersonic MHD turbulence 13 

Acknowledgments 

PG worked on section 3 and conducted the Enzo simulations. DV worked on section 
2. The analysis framework is a joint development by PG and DV. PG acknowledges 
financial support by the International Max Planck Research School for Solar System 
Science at the University of Gottingen. DV, WS and DS acknowledge research funding 
by the Deutsche Forschungsgemeinschaft (DFG) under grant SFB 963/1, project A15. 
CF thanks for funding provided by the ARC (grants DP130102078 and DP150104329). 
We gratefully acknowledge the Jlilich Supercomputing Centre (grant hhd20), the 
Leibniz Rechenzentrum and the Gauss Centre for Supercomputing (grant pr321o), the 
Partnership for Advanced Computing in Europe (PR ACE grant pr98mu), and the 
Australian National Computing Infrastructure (grant ek9). The software used in this 
work was in part developed by the DOE-supported Flash Center for Computational 
Science at the University of Chicago. 

References 

[1] Goldstein M L, Roberts D A and Matthaeus W H 1995 Ann. Rev. Astron. Astrophys. 33 283-326 

[2] Mac Low M M and Klessen R S 2004 Reviews of Modern Physics 76 125-194 

[3] Elmegreen B G and Scalo J 2004 Ann. Rev. Astron. Astrophys. 42 211-273 

[4] Cooper C M, Wallace J, Brookhart M, Clark M, Collins C, Ding W X, Flanagan K, Khalzov I, Li 

Y, Milhone J, Nornberg M, Nonn P, Weisberg D, Whyte D G, Zweibel E and Forest C B 2014 
Physics of Plasmas 21 013505 

[5] Schmidt W 2014 Living Reviews in Computational Astrophysics Submitted for publication 

[6] Gamier E, Adams N and Sagaut P 2009 Large Eddy Simulation for Compressible Flows 

[7] Muller W C and Carati D 2002 Physics of Plasmas (199f-present) 9 824-834 URL 

http://scitation.aip.org/content/aip/journal/pop/9/3/ 10.1063/1.1448498 

[8] Theobald M L, Fox P A and Sofia S 1994 Physics of Plasmas 1 3016-3032 

[9] Chernyshov A A, Karelsky I< V and Petrosyan A S 2014 Physics-Uspekhi 57 421-452 ISSN 1063- 

7869 URL http://stacks.iop.Org/1063-7869/57/5/421/ 

[10] Miki K and Menon S 2008 Phys. Plasmas 15 072306 ISSN 1070664X URL 

http://link.aip.org/link/PHPAEN/vl5/i7/p072306/sl&Agg=doi 

[11] Balarac G, Kosovichev A G, Brugiere O, Wray A A and Mansour N N 2010 Modeling of the 

subgrid-scale term of the filtered magnetic field transport equation Proceedings of the Summer 
Program (Center for Turbulence Research, Stanford University/NASA) pp 503-512 

[12] Favre A 1983 Physics of Fluids 26 2851 -2863 

[13] Chernyshov A A, Karelsky K V and Petrosyan A S 2007 Physics of Fluids (1994-present) 19 055106 

URL http://scitation.aip.Org/content/aip/journal/pof2/19/5/ 10.1063/1.2728936 

[14] Smagorinsky J 1963 Monthly Weather Review 91 99 

[15] Yoshizawa A 1990 Physics of Fluids B 2 1589-1600 

[16] Schmidt W and Federrath C 2011 Astron. Astrophys. 528 A106 ISSN 0004-6361 URL 

http://www.aanda.org/10.1051/0004-6361/201015630 

[17] Woodward P R, Porter D H, Anderson S, Fuchs T and Herwig F 2006 Journal of Physics 

Conference Series 46 370-384 

[18] Bryan G L, Norman M L, O’Shea B W, Abel T, Wise J H, Turk M J, Reynolds D R, Collins 

D C, Wang P, Skillman S W, Smith B, Harkness R P, Bordner J, hoon Kim J, Kuhlen M, Xu H, 
Goldbaum N, Hummels C, Kritsuk A G, Tasker E, Skory S, Simpson C M, Hahn O, Oishi J S, So 


Nonlinear closures for scale separation in supersonic MHD turbulence 


14 


G C, Zhao F, Cen R, Li Y and Collaboration T E 2014 The Astrophysical Journal Supplement 
Series 211 19 URL http://stacks.iop.org/0067-0049/211/i=2/a=19 

[19] Fryxell B, Olson K, Ricker P, Timmes F X, Zingale M, Lamb D Q, MacNeice P, Rosner R, 

Truran J W and Tufo H 2000 Astrophys. J. Suppl. Ser. 131 273-334 ISSN 0067-0049 URL 
http://stacks.iop.org/0067-0049/131/i=l/a=273 

[20] Waagan K, Federrath C and Klingenberg C 2011 Journal of Computational Physics 230 3331-3351 

[21] Federrath C, Chabrier G, Schober J, Banerjee R, Klessen R S and Schle¬ 

icher D R G 2011 Phys. Rev. Lett. 107 114504 ISSN 0031-9007 URL 

http://link.aps.org/doi/10.1103/PhysRevLett.107.114504 

[22] Federrath C, Schober J, Bovino S and Schleicher D R G 2014 Astrophys. J. 

797 L19 ISSN 2041-8213 URL http://stacks. iop. org/2041-8205/797/i=2/ 

a=L19?key=crossref.2bfd2aece4bf89209f7a97219f506f70 

[23] Dedner A, Kemm F, Kroner D, Munz C D, Schnitzer T and Wesenberg M 

2002 Journal of Computational Physics 175 645 - 673 ISSN 0021-9991 URL 

http://www.sciencedirect.com/science/article/pii/ S00219991019696IX 

[24] Schmidt W, Federrath C, Hupp M, Kern S and Niemeyer J C 2009 Astronomy & Astrophysics 

494 127-145 

[25] Federrath C, Roman-Duval J, Klessen R S, Schmidt W and Mac Low M M 2010 

Astron. & Astrophys. 512 A81 

[26] Federrath C 2013 Mon. Not. R. Astron. Soc. 436 1245-1257 

ISSN 0035-8711 URL http://arxiv.org/abs/1306.3989 http:// 

mnras.oxfordj ournals.org/cgi/doi/10.1093/mnras/stt1644 

[27] Galtier S and Banerjee S 2011 Phys. Rev. Lett. 107 134501 ISSN 0031-9007 URL 

http://link.aps.org/doi/10.1103/PhysRevLett.107.134501 

[28] Banerjee S and Galtier S 2013 Phys. Rev. E 1-8 URL 

http://pre.aps.org/abstract/PRE/v87/il/e013019 http:// arxiv.org/abs/1301.2470 

[29] Aluie H 2011 Phys. Rev. Lett. 106 174502 ISSN 0031-9007 URL 

http://link.aps.org/doi/10.1103/PhysRevLett.106.174502 

[30] Aluie H 2013 Physica D Nonlinear Phenomena 247 54-65 

[31] Kitsionas S, Federrath C, Klessen R S, Schmidt W, Price D J, Dursi L J, Gritschneder M, Walch 

S, Piontek R, Kim J, Jappsen A K, Ciecielag P and Mac Low M M 2009 Astron. Astrophys. 
508 541-560 ISSN 0004-6361 URL http://www.aanda.org/10.1051/0004-6361/200811170 

[32] Clark P C, Glover SCO, Smith R J, Greif T H, Klessen R S and Bromm V 2011 Science 331 

1040- 

[33] Greif T H, Schleicher D R G, Johnson J L, a K Jappsen, Klessen R S, Clark P C, Glover 

SCO, Stacy a and Bromm V 2008 Proc. Int. Astron. Union 4 33 ISSN 1743-9213 URL 
http://www.journals.Cambridge.org/ abstract_S1743921308024551 

[34] Greif T H, Bromm V, Clark P C, Glover SCO, Smith R J, Klessen R S, Yoshida N and 

Springel V 2012 Formation and evolution of primordial protostellar systems vol 51 pp 51- 
56 ISBN 9780735410923 URL http://scitation.aip.org/content/aip/proceeding/aipcp/ 
10.1063/1.4754327 

[35] Agertz O, Kravtsov A V, Leitner S N and Gnedin N Y 2013 Astrophys. J. 

770 25 ISSN 0004-637X URL http://stacks. iop. org/0004-637X/770/i=l/ 

a=25?key=crossref.20941ed3bll63ecbe48734ff39a64be2 

[36] Devriendt J, Rimes C, Pichon C, Teyssier R, Le Borgne D, Aubert D, Audit E, Colombi S, Courty 

S, Dubois Y, Prunet S, Rasera Y, Slyz A and Tweed D 2010 Mon. Not. R. Astron. Soc. Lett. 
403 L84-L88 ISSN 17453925 URL http://mnrasl.oxfordjournals.org/cgi/doi/10.llll/ 
j.1745-3933.2010.00820.x 

[37] Vazza F, Briiggen M, Gheller C and Wang P 2014 000 URL http: //arxiv. org/abs/1409.2640v2 

http://arxiv.org/abs/ 1409.2640 

[38] Schleicher D, Latif M, Schober J, Schmidt W, Bovino S, Federrath C, Niemeyer J, 


Nonlinear closures for scale separation in supersonic MHD turbulence 


15 


Banerjee R and Klessen R 2013 Astron. Nachrichten 334 531-536 ISSN 00046337 URL 
http://doi.wiley.com/10.1002/asna.201211898 

[39] Latif M A, Schleicher D R G, Schmidt W and Niemeyer J 2013 Mon. Not. R. Astron. Soc. 432 

668-678 ISSN 0035-8711 URL http: //mnras. oxford journals. org/cgi/doi/10.1093/mnras/ 
stt503 

[40] Latif M A, Schleicher D R G and Schmidt W 2014 Mon. Not. R. Astron. Soc. 440 1551 -1561 ISSN 

0035-8711 URL http://mnras.oxfordjournals.org/cgi/doi/10.1093/mnras/ stu357 


